#include "../../Research2.0/Modeling/GernAlgSurf.cuh"


__global__ void KIOneSolve(unsigned char* Normal, unsigned char* Pos, size_t pitch, uint degSum, GAS_SURF gas)
{
	uint x = blockIdx.x*blockDim.x + threadIdx.x;
	uint y = blockIdx.y*blockDim.y + threadIdx.y;
	if(x>=cRenderParam.renderWidth || y>=cRenderParam.renderHeight)			return;

	float* pixel = (float*)(Pos + y*pitch) + 4*x;
	pixel[0] = 1.0;
	pixel[1] = 1.0;
	pixel[2] = 1.0;
	pixel[3] = 0.0f;


	float dir[3] ={2*(x+0.0f)/cRenderParam.renderWidth-1.0f, 1.0f-2*(y+0.0f)/cRenderParam.renderHeight, 0.0f};

	float trans[4];
	trans[0] = dir[0]*cRenderParam.invViewProj[0][0]+dir[1]*cRenderParam.invViewProj[1][0]+dir[2]*cRenderParam.invViewProj[2][0]+cRenderParam.invViewProj[3][0]; 
	trans[1] = dir[0]*cRenderParam.invViewProj[0][1]+dir[1]*cRenderParam.invViewProj[1][1]+dir[2]*cRenderParam.invViewProj[2][1]+cRenderParam.invViewProj[3][1]; 
	trans[2] = dir[0]*cRenderParam.invViewProj[0][2]+dir[1]*cRenderParam.invViewProj[1][2]+dir[2]*cRenderParam.invViewProj[2][2]+cRenderParam.invViewProj[3][2]; 
	trans[3] = dir[0]*cRenderParam.invViewProj[0][3]+dir[1]*cRenderParam.invViewProj[1][3]+dir[2]*cRenderParam.invViewProj[2][3]+cRenderParam.invViewProj[3][3]; 
	dir[0] = trans[0]/trans[3]-cRenderParam.eye[0];
	dir[1] = trans[1]/trans[3]-cRenderParam.eye[1];
	dir[2] = trans[2]/trans[3]-cRenderParam.eye[2];

	float tmin, tmax;
	int firstSeg[3] = {0, 0, 0};
	int seg[3] = {1, 1, 1};
	if(0==RayBoxIntersect(&tmin, &tmax, cRenderParam.eye, dir, gas.first, gas.def, firstSeg, seg)|| tmax<1 ||fabsf(tmin)>1.0e5||fabsf(tmax)>1.0e5)	return;
	if(tmin<1)	tmin=1;


	float fValue[MAX_DEG_INTERPOLATION+1];
#if _KI_BEZ_
	float coeff[MAX_DEG_INTERPOLATION+1];
#else
	float coeff[MAX_CONTROL_POINT];
#endif

	//********* code surround by /**/ is chebyshev interpolation ****************//
	//remember: change the interpolation here will also need change the InvMatrix;
	float In[3] ={cRenderParam.eye[0]+tmin*dir[0],cRenderParam.eye[1]+tmin*dir[1],cRenderParam.eye[2]+tmin*dir[2]};
	float Out[3]={cRenderParam.eye[0]+tmax*dir[0],cRenderParam.eye[1]+tmax*dir[1],cRenderParam.eye[2]+tmax*dir[2]};
	float xyz[3];
	float cheby;
	for(int i=0; i<=degSum; i++)
	{
		cheby = (1-(__cosf(PI*(2*i+1)/(2*(degSum+1)))))/2;
		xyz[0] = (1-cheby)*In[0]+cheby*Out[0];
		xyz[1] = (1-cheby)*In[1]+cheby*Out[1];
		xyz[2] = (1-cheby)*In[2]+cheby*Out[2];
		fValue[i] = EvalGenAlgSurf(xyz[0], xyz[1], xyz[2]);
	}
	//********* code surround by /**/ is chebyshev interpolation ****************//


	//********* code surround by /**/ is equidistance interpolation ****************//
	//float tstep = (tmax-tmin)/degSum;
	//float xyz[3]={(cRenderParam.eye[0]+tmin*dir[0]), (cRenderParam.eye[1]+tmin*dir[1]), (cRenderParam.eye[2]+tmin*dir[2])};
	//dir[0] *= tstep;
	//dir[1] *= tstep;
	//dir[2] *= tstep;
	//for(int i=0; i<=degSum; i++)
	//{
	//	fValue[i] = EvalGenAlgSurf(xyz[0], xyz[1], xyz[2]);
	//	xyz[0] += dir[0];
	//	xyz[1] += dir[1];
	//	xyz[2] += dir[2];
	//}
	//dir[0] /= tstep;
	//dir[1] /= tstep;
	//dir[2] /= tstep;
	//********* code surround by /**/ is equidistance interpolation ****************//


	uint matrixInd = 0;
	for(int i=0; i<=degSum; i++)
	{
		coeff[i] = 0.0;
		for(int j=0; j<=degSum; j++)
		{
			coeff[i] += fValue[j]*cInvMatrix[matrixInd++];
		}
	}
	int indC;
	for(indC=0; indC<degSum; indC++)
	{
		if(coeff[indC]*coeff[indC+1]<0)	break;
	}
	if(indC==degSum)	return;

	if(cRenderParam.material.transparent)
	{	
		float root[MAX_DEG_INTERPOLATION];
#if _KI_BEZ_
		int rootN = KnotInsertBezAll(coeff, degSum, root);
#else
		int rootN = KnotInsertAll(coeff, degSum, root);
#endif
		if(rootN<=0)			return;

		float rt;
		float alphaB=cRenderParam.material.t, alphaBB=1.0f-alphaB;
		float* pixel = (float*)(Pos+y*pitch)+4*x;
		float rgb[3] = {0,0,0}; 
		float matDiff[3], Pos[3], ray[3], diffuse, norm[3];
		matDiff[0] = cRenderParam.material.d[0];
		matDiff[1] = cRenderParam.material.d[1];
		matDiff[2] = cRenderParam.material.d[2];	

		for(int i=rootN-1;i>=0;i--)
		{
			rt = root[i];
			if(rt<=0||rt>=1.0)	continue;
			rt = tmin*(1-rt) + tmax*rt;	

			Pos[0] = cRenderParam.eye[0]+rt*dir[0];
			Pos[1] = cRenderParam.eye[1]+rt*dir[1];
			Pos[2] = cRenderParam.eye[2]+rt*dir[2];
			EvalGenAlgSurfDf(Pos[0], Pos[1], Pos[2], norm);
			normalize(norm);


			ray[0] = cRenderParam.eye[0]-Pos[0];
			ray[1] = cRenderParam.eye[1]-Pos[1];
			ray[2] = cRenderParam.eye[2]-Pos[2];
			normalize(ray);
			diffuse = fabs(norm[0]*ray[0]+norm[1]*ray[1]+norm[2]*ray[2]);
			diffuse *= diffuse;
			//diffuse *= diffuse;
			diffuse *= 3;


			rgb[0] = alphaBB*(diffuse*matDiff[0]) + alphaB*rgb[0];
			rgb[1] = alphaBB*(diffuse*matDiff[1]) + alphaB*rgb[1];
			rgb[2] = alphaBB*(diffuse*matDiff[2]) + alphaB*rgb[2];
		}
		pixel[0] = rgb[0];
		pixel[1] = rgb[1];
		pixel[2] = rgb[2];
		pixel[3] = 1.0f;
	}
	else
	{
 
 #if _KI_BEZ_
  		float firstRoot = KnotInsertBez(coeff, degSum);
 #else
  		float firstRoot = KnotInsert(coeff, degSum);
 #endif
 
  		if(firstRoot<0)		return;
 		firstRoot = tmin*(1-firstRoot) + tmax*firstRoot;
 		xyz[0] = (cRenderParam.eye[0]+firstRoot*dir[0]);
 		xyz[1] = (cRenderParam.eye[1]+firstRoot*dir[1]);
 		xyz[2] = (cRenderParam.eye[2]+firstRoot*dir[2]);
 		pixel[0] = xyz[0];
 		pixel[1] = xyz[1];
 		pixel[2] = xyz[2];
 		pixel[3] = firstRoot;


		//for polar surface demo
//		float root[MAX_DEG_INTERPOLATION];
//#if _KI_BEZ_
//		int rootN = KnotInsertBezAll(coeff, degSum, root);
//#else
//		int rootN = KnotInsertAll(coeff, degSum, root);
//#endif
//		if(rootN==0)	return;
//		float rt;
//		xyz[0] = xyz[1] = xyz[2] = 100;
//		for(int i=0; i<rootN; i++)
//		{
//			rt = root[i];
//			if(rt<0||rt>1.0)	break;
//			rt = tmin*(1-rt) + tmax*rt;
//			xyz[0] = (cRenderParam.eye[0]+rt*dir[0]);
//			xyz[1] = (cRenderParam.eye[1]+rt*dir[1]);
//			xyz[2] = (cRenderParam.eye[2]+rt*dir[2]);
//			if(xyz[2]<(cRenderParam.nIter-5)/20.0)	break;
//		}
//		if(xyz[2]>(cRenderParam.nIter-5)/20.0||xyz[2]==100)	
//		{
//			pixel[0] = pixel[1] = pixel[2] = pixel[3] = 0;
//		}
//		else
//		{
//			pixel[0] = xyz[0];
//			pixel[1] = xyz[1];
//			pixel[2] = xyz[2];
//			pixel[3] = rt;
//		}



#if	COMPUTE_NORMAL
		pixel = (float*)(Normal + y*pitch) + 4*x;
		float norm[3];
		EvalGenAlgSurfDf(xyz[0], xyz[1], xyz[2], norm);
		//norm[0] = EvalGenAlgSurf(xyz[0]+DERIVATIVE_DELTA, xyz[1], xyz[2]);
		//norm[1] = EvalGenAlgSurf(xyz[0], xyz[1]+DERIVATIVE_DELTA, xyz[2]);
		//norm[2] = EvalGenAlgSurf(xyz[0], xyz[1], xyz[2]+DERIVATIVE_DELTA);
		//norm[0] *= DERIVATIVE_DELTA_INV;
		//norm[1] *= DERIVATIVE_DELTA_INV;
		//norm[2] *= DERIVATIVE_DELTA_INV;

		normalize(norm);
		pixel[0] = norm[0];
		pixel[1] = norm[1];
		pixel[2] = norm[2];
#endif
	}

	//float rgb[3] = {0.0, 0.0, 0.0};
	////MATERIAL material;
	////chooseMaterial(0, &material, &cRenderParam.material);
	////woodMaterial(xyz, &material, &cRenderParam.material);
	//shading(rgb, xyz, norm, cRenderParam.eye, cRenderParam.nLight, cRenderParam.light, &cRenderParam.material, 0);
	//pixel[0] = rgb[0];
	//pixel[1] = rgb[1];
	//pixel[2] = rgb[2];
}


extern "C" void KIOneKer(RENDER_PARAM* renderParam, unsigned char* surfNormal, unsigned char* surface, size_t pitch, bool* reNormal)
{
	cudaMemcpyToSymbol(cRenderParam, renderParam, sizeof(RENDER_PARAM));
	
	CUDA_SAFE_CALL( cudaMemset(surfNormal, 0, pitch*height) );
	CUDA_SAFE_CALL( cudaMemset(surface, 0, pitch*height) );

	dim3 Db = dim3(8,8);
	dim3 Dg = dim3( width/Db.x+1, height/Db.y+1 );
	if(width %Db.x==0)	Dg.x-=1;
	if(height%Db.y==0)  Dg.y-=1;

#if	COMPUTE_NORMAL
	*reNormal = true;
#else
	//if(renderParam->material.transparent==1||renderParam->renderState==1)
	//	*reNormal = true;
	//else
		*reNormal = false;
#endif

	KIOneSolve<<<Dg, Db>>>((unsigned char*)surfNormal, (unsigned char*)surface, pitch, degSum, GASData);
}

